clear all
close all
clc
%Change rxn rate on line 43 to go from instantaneous 1e-5 to very slow 1e-10

%% D/H Isotope Diffusion-Reaction Model
% for the glass-vapor system

beta = 1; %
kf = 1E-7; %1e-4 = fast reaction; 1e-10 slow reaction
alpha_oh_m = 0.96; %fractionation between OH-H2Om

%D/H diffusion - reaction model
TC = 225; %T in Celsius
TK = TC + 273.15; %T in Kelvin
total_time = 500; %time in hours
totalt = total_time*60*60; %hours converted to seconds

Keq = exp(1.89-3120/TK);
kb = kf/Keq;
alpha_oh_m = alpha_oh_m;%alpha_OH/H2Om; reaction kinetics
alpha_f = 0.5;%kinetic fractionation factor for H2Om + O --> 2OH
af = alpha_f*kf;
ab = af/(Keq*alpha_oh_m);
beta = beta;%Prefactor for HDOm diffusion coefficient (1 = same as DH2Om)

vsmow = 1/6420;
delta_vapor = 75.3; %permil
delta_gb = 50.3;%Composition of total water at the boundary before reaction; 
%for instantaneous reaction this value remains at +50. for slow reaction,...
%the composition of H2Ot matched H2Om which is heavier than +50.
delta_glass_initial = -100.8;
Rti = (delta_glass_initial/1000 + 1)*vsmow;
% P is T-dependent in these experiments
% 175C = 0.89MPa; 225C = 2.54MPa; 275C = 5.92MPa; 375C = 21.1MPa
P = 2.54; %P in MPa
pf = 2; % DH2O prefactor (Hudak + Bindeman, 2020)

lengthx = 160;% length of domain in microns
xnodes = 200;
dx = lengthx/(xnodes-1);
nodes = zeros(xnodes,1);
for i = 2:xnodes
    nodes(i) = nodes(i-1)+dx;
end
for i = 1:xnodes-1
    Vshell(i) = 4/3*pi*(nodes(i+1)^3) - 4/3*pi*(nodes(i)^3);
end
volume = sum(Vshell);
%%
%Arrays
H2Ot = zeros(xnodes,1); %H2Ot in mole fraction
wH2Ot = zeros(xnodes,1); %H2Ot in wt%
XH2Om = zeros(xnodes,1);
H2Om = zeros(xnodes,1);
H2Om_new = zeros(xnodes,1);
wH2Om = zeros(xnodes,1);
OH = zeros(xnodes,1);
wOH = zeros(xnodes,1);
O = zeros(xnodes,1);
OH_new = zeros(xnodes,1);

HDOm = zeros(xnodes,1);
wHDOm = zeros(xnodes,1);
HDOm_new = zeros(xnodes,1);
OD = zeros(xnodes,1);
wOD = zeros(xnodes,1);
OD_new = zeros(xnodes,1);

dOH = zeros(xnodes,1);
dOD = zeros(xnodes,1);
dH2Om = zeros(xnodes,1);
dHDOm = zeros(xnodes,1);

d2OH = zeros(xnodes,1);
d2OD = zeros(xnodes,1);
d2H2Om = zeros(xnodes,1);
d2HDOm = zeros(xnodes,1);

DH2Om = zeros(xnodes,1); %um^2/s
DHDOm = zeros(xnodes,1);
dD_H2Om = zeros(xnodes,1);
dD_HDOm = zeros(xnodes,1);

Rt = zeros(xnodes,1);
ROH = zeros(xnodes,1);
RH2Om = zeros(xnodes,1);
deltat = zeros(xnodes,1);
delta_OH = zeros(xnodes,1);
delta_H2Om = zeros(xnodes,1);

step = 0;
%%
%This section has been validated against H2O_speciation excel spreadsheet
wH2Oti = 0.1;
Cti = wH2Oti/100;
H2Oti = (Cti/18.015)/(Cti/18.015+(1-Cti)/32.49);%eq. 2.81 Zhang 2008
H2Omi = 8*H2Oti^2/(Keq*(1-2*H2Oti)+8*H2Oti+(Keq^2*(1-2*H2Oti)^2+16*Keq*H2Oti*(1-H2Oti))^0.5);%eq. 3.83a Zhang 2008
OHi = 2*(H2Oti-H2Omi);%eq. 2.83 Zhang 2008
Oi = 1-H2Omi-OHi; %2.84 %TREAT OD AND D2O AS TRACER (LIKE CHEN MODEL)
a = (0.5/alpha_oh_m)*af/ab*Oi/(OHi^2);
b = 0.5*af*Oi/(ab*OHi)+1-Rti;
c = -Rti*(2*H2Omi+OHi);
p = [a b c];
r1 = roots(p);
HDOmi = max(r1);
%HDOmi = Rti*(2*H2Omi+OHi)/(1+alpha_eq/2*OHi/H2Omi-Rti);%APPROXIMATION
ODi = alpha_oh_m*(0.5*HDOmi/H2Omi)*OHi;

% H2O solubility in glass = 3-5wt% for 175-275C (Hudak + Bindeman, 2020)
wH2Otb = 3.9;
Ctb = wH2Otb/100;
Rtb = (delta_gb/1000 + 1) * vsmow;
H2Otb = (Ctb/18.015)/(Ctb/18.015+(1-Ctb)/32.49);%eq. 2.81, Zhang 2008
H2Omb = 8*H2Otb^2/(Keq*(1-2*H2Otb)+8*H2Otb+(Keq^2*(1-2*H2Otb)^2+16*Keq*H2Otb*(1-H2Otb))^0.5);%eq. 3.83a Zhang 2008
OHb = 2*(H2Otb-H2Omb);%eq. 2.83 Zhang 2008
Ob = 1-H2Omb-OHb; %2.84
a = (0.5/alpha_oh_m)*af/ab*Ob/(OHb^2);
b = 0.5*af*Ob/(ab*OHb)+1-Rtb;
c = -Rtb*(2*H2Omb+OHb);
p = [a b c];
r2 = roots(p);
HDOmb = max(r2);
%HDOmb = Rtb*(2*H2Omb+OHb)/(1+alpha_eq/2*OHb/H2Omb-Rtb);%APPROXIMATION
ODb = alpha_oh_m*(0.5*HDOmb/H2Omb)*OHb;
%%
for i = 1:xnodes
    H2Ot(i) = H2Oti;
    H2Om(i) = H2Omi;
    OH(i) = OHi;
    O(i) = Oi;
    HDOm(i) = HDOmi;
    OD(i) = ODi;
    Rt(i) = Rti;
end
%%
%Perturbation at glass-vapor interface
H2Om(end) = H2Omb;
HDOm(end) = HDOmb;
O(end) = Ob;
Rt(end) = Rtb;
%Calculate the appropriate time step from the max diffusivity and create arrays to store the cumulative H2O and RT of the glass
DH2Ommax = exp((14.08-13128/TK-2.796*P/TK)+H2Otb*(-27.21+36892/TK+57.23*P/TK))*pf; %*(H2Om(i)/(H2Om(i)+HDOm(i)));
t = 0; %seconds
dt = (0.1*dx^2)/DH2Ommax;%
nsteps = ceil(totalt/dt);
time = zeros(nsteps,1);
H2Obulk = zeros(nsteps,1);
dDbulk = zeros(nsteps,1);
alpha_bulk = zeros(nsteps,1);
alpha_OH = zeros(nsteps,1);
alpha_H2Om = zeros(nsteps,1);
deltat_end = zeros(nsteps,1);
deltaOH_end = zeros(nsteps,1);
deltaH2Om_end = zeros(nsteps,1);

%%               1 D   S P H E R I C A L   D I F F U S I O N
step=1;
while t < totalt
    %Update diffusivities
    for i = 1:xnodes
        DH2Om(i) = exp((14.08-13128/TK-2.796*P/TK)+H2Ot(i)*(-27.21+36892/TK+57.23*P/TK))*pf;
        DHDOm(i) = beta*DH2Om(i);
    end
    %Spatial derivatives
    for i = 2:xnodes-1
        dH2Om(i) = (H2Om(i+1)-H2Om(i-1))/(2*dx);
        dHDOm(i) = (HDOm(i+1)-HDOm(i-1))/(2*dx);
        dOH(i) = (OH(i+1)-OH(i-1))/(2*dx);
        dOD(i) = (OD(i+1)-OD(i-1))/(2*dx);
        dD_H2Om(i) = (DH2Om(i+1)-DH2Om(i-1))/(2*dx);
        dD_HDOm(i) = (DHDOm(i+1)-DHDOm(i-1))/(2*dx);
        d2H2Om(i) = (H2Om(i+1)-2*H2Om(i)+H2Om(i-1))/(dx^2);
        d2HDOm(i) = (HDOm(i+1)-2*HDOm(i)+HDOm(i-1))/(dx^2);
        d2OH(i) = (OH(i+1)-2*OH(i)+OH(i-1))/(dx^2);
        d2OD(i) = (OD(i+1)-2*OD(i)+OD(i-1))/(dx^2);
    end
    %Diffusion
    for i = 2:xnodes-1
        %        H2Om(i) = (dD_H2Om(i)*dH2Om(i)+DH2Om(i)*d2H2Om(i))*dt + H2Om(i);
        %        HDOm(i) = (dD_HDOm(i)*dHDOm(i)+DHDOm(i)*d2HDOm(i))*dt + HDOm(i);
        H2Om(i) = (2/nodes(i)*DH2Om(i)*dH2Om(i)+dD_H2Om(i)*dH2Om(i)+DH2Om(i)*d2H2Om(i))*dt + H2Om(i);
        HDOm(i) = (2/nodes(i)*DHDOm(i)*dHDOm(i)+dD_HDOm(i)*dHDOm(i)+DHDOm(i)*d2HDOm(i))*dt + HDOm(i);
        %        OH(i) = (2/nodes(i)*DH2Om(i)*dOH(i)+dD_H2Om(i)*dOH(i)+DH2Om(i)*dOH(i))+OH(i);
        %        OD(i) = (2/nodes(i)*DHDOm(i)*dOD(i)+dD_HDOm(i)*dOD(i)+DHDOm(i)*dOD(i))+OD(i);
    end
    %Reaction
    for i=1:xnodes
        H2Om_new(i) = (-kf*H2Om(i)*O(i)+kb*OH(i)^2)*dt + H2Om(i);
        OH_new(i) = (kf*H2Om(i)*O(i)-kb*OH(i)^2)*dt + OH(i);
        HDOm_new(i) = 2*(-0.5*af*HDOm(i)*O(i)+ab*OD(i)*OH(i))*dt + HDOm(i);
        OD_new(i) = (0.5*af*HDOm(i)*O(i)-ab*OD(i)*OH(i))*dt + OD(i);
        %O(i) = (-kf*H2Om(i)*O(i)+kb*OH(i)^2)*dt+O(i);
    end
    for i = 1:xnodes
        H2Om(i) = H2Om_new(i);
        HDOm(i) = HDOm_new(i);
        OH(i) = OH_new(i);
        OD(i) = OD_new(i);
    end
    %No-flux BC at center of particle
    H2Om(1) = H2Om(2);
    OH(1) = OH(2);
    O(1) = O(2);
    HDOm(1) = HDOm(2);
    OD(1) = OD(2);
    %Update H2Ot (for diffusivity) and delta values (for plotting)
    for i = 1:xnodes
        H2Ot(i) = H2Om(i)+OH(i)/2; %mole fraction on a single oxygen basis
        C(i) = (-H2Ot(i)/32.49)/((H2Ot(i)/18.015-H2Ot(i)/32.49-1/18.015));
        wH2Ot(i) = 100*C(i);%wt%
        XH2Om(i) = H2Om(i)/(H2Om(i)+OH(i));
        O(i) = 1-H2Om(i)-OH(i);
        delta_OH(i) = 1000*(OD(i)/OH(i)/vsmow - 1);
        delta_H2Om(i) = 1000*(0.5*HDOm(i)/H2Om(i)/vsmow - 1);
        %Rt(i) = XH2Om(i)*RH2Om(i)+(1-XH2Om(i))*ROH(i);
        Rt(i) = (HDOm(i)+OD(i))/(2*H2Om(i)+HDOm(i)+OH(i));%APPROXIMATION
        deltat(i) = 1000*(Rt(i)/vsmow - 1);
    end
    H2Om(end) = H2Omb;
    HDOm(end) = HDOmb;

    %Calculate and store the bulk H2O and RT of the glass
    for i = 1:xnodes-1
        H2O_weighted(i) = (wH2Ot(i)+wH2Ot(i+1))/2 * Vshell(i);
        dD_weighted(i) = (wH2Ot(i)*deltat(i) + wH2Ot(i+1)*deltat(i+1))/2 * Vshell(i);
    end
    H2Obulk(step) = sum(H2O_weighted)/volume;
    dDbulk(step) = sum(dD_weighted)/(volume*H2Obulk(step));
    time(step) = t/3600;
    alpha_bulk(step) = (1000+deltat(end))/(1000+delta_vapor);
    alpha_OH(step) = (1000+delta_OH(end))/(1000+delta_H2Om(end));
    alpha_H2Om(step) = (1000+delta_H2Om(end))/(1000+delta_vapor);
    deltat_end(step) = deltat(end);
    deltaOH_end(step) = delta_OH(end);
    deltaH2Om_end(step) = delta_H2Om(end);
    t = t + dt;
    step = step +1;
end
%%
%Convert from mol fraction on a single oxygen basis to wt%
for i = 1:xnodes
    C(i) = (-H2Ot(i)/32.49)/((H2Ot(i)/18.015-H2Ot(i)/32.49-1/18.015));
    wH2Ot(i) = 100*C(i);
    wH2Om(i) = 100*(H2Om(i)*C(i)/H2Ot(i)); %eq. 2.82 (wt%)
    wOH(i) = wH2Ot(i)-wH2Om(i);%top of p. 125 (wt%)
end

%% figures

figure(1)
subplot(1,2,1)
plot(time,H2Obulk,'-')
xlabel('Time (hours)')
ylabel('H_2O (wt.%)')
xlim([0 total_time])

subplot(1,2,2)
plot(time,dDbulk,'-')
xlabel('Time (hours)')
ylabel('\deltaD (permil)')
ylim([-101 50])
xlim([0 total_time])

